Simulation of wavepacket tunneling of interacting identical particles 
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We demonstrate a new method of simulation of nonstationary quantum processes, considering the 
tunneling of two interacting identical particles, represented by wave packets. The used method of 
quantum molecular dynamics (WMD) is based on the Wigner representation of quantum mechanics. 
In the context of this method ensembles of classical trajectories are used to solve quantum Wigner- 
Liouville equation. These classical trajectories obey Hamilton-like equations, where the effective 
potential consists of the usual classical term and the quantum term, which depends on the Wigner 
function and its derivatives. The quantum term is calculated using local distribution of trajectories 
in phase space, therefore classical trajectories are not independent, contrary to classical molecular 
dynamics. The developed WMD method takes into account the influence of exchange and interaction 
between particles. The role of direct and exchange interactions in tunneling is analyzed. The 
tunneling times for interacting particles are calculated. 

PACS numbers: 34.10.+X, 03.65.Xp, 71.15.Pd, 02.70.Ns 



I. INTRODUCTION 



A new quantum molecular dynamics method (QMD) 
was recently successfully applied to a single wavepacket 
tunneling [l|, This method is based on the Wigner 
representation |3, Q of quantum mechanics (further ref- 
ered to as WMD - the Wigner representation based MD) . 
In the present paper we further develop this method 
and consider its application to the many-body problem of 
nonstationary tunneling of interacting identical particles. 
Nonstationary tunneling is a problem of great interest in 
particular in connection with developments of nanoelec- 
tronics. Until now role of interaction and exchange in 
nonstationary tunneling is not clear. To clear up this 
question is one the aims of this paper. In this connec- 
tion we consider the tunneling of two identical charged 
particles, represented by wavepackets. 

In the Wigner representation of quantum mechanics 
the state of the system is described by the Wigner func- 
tion, which obeys Wigner-Liouville equation. The equa- 
tion can be rewritten in the form analogous to classi- 
cal Liouville equation for classical distribution function. 
This analogy is the basis of WMD: the ensembles of clas- 
sical trajectories are used to solve numerically quantum 
Wigner-Liouville equation. The trajectories can be de- 
termined by equations of motion analogous to classical 
ones. The used modification as against classical equa- 
tions of motion for the trajectories is an addition of extra 
'quantum' term in the expression for the force 0. This 
'quantum' term is expressed through the local approxi- 
mation of the Wigner function. For the approximation 
of the Wigner function we used multi-dimensional Gauss 



distribution with the parameters determined through the 
local moments of the ensemble of classical trajectories. 

In the present paper the wavepackets moving in 
double-well potential were considered. The interparticle 
interactions are fully taken into account. The wavepack- 
ets are initially placed in the same well on the one side 
of the barrier. We analyzed the long-time evolution of 
wavepackets (for time scales corresponding to many pe- 
riods oscillation in the well) and consider the probability 
to detect a particle in the first and in the second well, 
respectively. Besides we study the short-time evolution 
(characteristic times of interaction of wavepacket with 
the barrier) and regard tunneling times. 

Tunneling time is one of the most important features 
of nonstationary tunneling. However, the theoretical def- 
inition of this quantity is nontrivial. There are exist a lot 
of definitions of tunneling time 

J,.3„ .14 .15, J£] . We use two common approaches to deter- 
mine tunneling time, namely presence and arrival times 
(see 0, and references therein) . 

First, one can consider the detector which reacts to the 
presence of particles at some point xq. The values mea- 
sured by this detector in a set of experiments on, e.g., 
particles transmission through a barrier, would depend 
on particle density p{xQ,t) at time t and the mean pres- 
ence time of a particle at point xq would be given by 



(tpi^o)) ^ / dtt p{xi),t) j / dtp{xo,t) (1) 



For two points Xa and Xb one can consider the average 
time of transmission: 



{tT{Xa,Xb)) = {tp{Xb)) ~ (tpiXa)) 



(2) 
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If the points are located on the different sides of a bar- 
rier, then expression |2Jl is an approximation for tunnel- 
ing time. 
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Second, the detector that measures flux density can 
be used. For this set of experiments one need to define 
another quantity - arrival time. In this case the flux 
density operator must be considered 



J{xq) ^[P^iQ - ^o) + S{q - Xo)p] 



m)\JM\m) 
dt{m\JM\m) 



J{xQ,t) 



and the arrival time at point a;o can be defined as 



(3) 
(4) 



dynamics (MD) , the distinction is only in the calculation 
of the force and in the probability interpretation of ini- 
tial conditions. During about 40 years the classical MD 
methods were sufficiently improved and all advantageous 
numerical schemes can be simply implemented in WMD. 
The modern MD techniques allows to operate with the 
thousands of particles and the same can be in principle 
achieved by means of WMD, but in the last case one can 
consider quantum particles. 

We describe the simulation method and the physical 
model in SecE] and Sec lIIII respectively. Results are 
presented and discussed in Sec II VI Main conclusions are 
summarized in SecfVl 



{ta{xo))^ I dttJ{xo,t)/ I dtJ{xo,t) 



(5) 



We stress that J{xQ,t) can be negative due to the op- 
posite flux. Therefore Eq. can be used directly as a 
probability distribution of arrival times only if the op- 
posite flux is negligible. This requirement can be ful- 
flUcd if the detector is located far from the barrier. Then 
Eq. Q determines the quasidistribution of the arrival 
times. But in this case one can not distinguish the time 
of transmission under the barrier and time of passing the 
region between the barrier and the detector - still un- 
resolved problem of time measurement in quantum me- 
chanics (see, e.g., 11]). We use the presence and arrival 
times from all variety of possible definitions of tunnel- 
ing time because their measurement in the framework of 
WMD is relatively simple and, what is more important, 
the physical meaning of the Eq. |^ and 10} is transpar- 
ent and connected with the use of point-like detectors in 
the set of the experiments on particles transmission. 

By changing the strength of interaction between the 
particles, we investigate the influence of interaction on 
tunneling. We also consider the role of exchange. We 
found that the exchange is important if the interaction is 
weak. In this case exchange has a substantial influence 
on both the tunneling probability and tunneling time. 
With the increase of interaction initial system energy 
with flxed initial wave functions becomes greater. This 
leads to decrease of tunneling times, the role of exchange 
gets smaller and tunneling becomes insignificant in com- 
parison with the passing above the barrier. Our inves- 
tigation had shown that WMD is advantageous method, 
which can be used to solve the many-body problems with- 
out enormous computer resources, and which allows to 
take into account such essentially quantum features as 
exchange and tunneling. 

We present here the investigation of the two-particle 
problem, but the generalization of WMD for the case of 
more particles is straightforward. The advantage of using 
the Wigner representation in comparison with direct nu- 
merical solution of Schrodinger equation is the following. 
Using WMD one does not need to store large number 
arrays as with the grid methods. The basic algorithm 
of WMD is very close to that of the common molecular 



II. SIMULATION METHOD 

A. Equations of motion for Wigner trajectories 

The Wigner representation of quantum mechanics is 
one of the representations which uses quantum distri- 
bution function in phase space. The Wigner function 
F^{q,p,t) describes time evolution of the system and 
average values of physical quantities are calculated with 
the help of scalar functions, Weyl symbols {q,p): 

{A) ^JdpJ dqA'^{q,p)F'^{q,p,t) (6) 

It can be shown 0, , that Weyl symbols are expressed 
through corresponding operators A{q,p) as follows: 



h 

2^ 



A^{q,p)^^ I dCd7]Tr[A{q,p)e'^^+"'P] g-'^'-^^P 



(7) 

The Wigner function is real and satisfies the following 
rules 



dpF'^{q,p,t) = {q\p\q), 



dqF'^{q,p,t)^{p\p\p), 



(8) 
(9) 



here p is the density operator. The Wigner function 
F^ {q,p,t) is also not nonnegative. There are nonnega- 
tive quantum distribution functions, for example Husimi 
function |l9j| . but its evolution equation is usually more 
complicated as against the Wigner function. 

If one considers the Hamiltonian H = p^ / (2m) + V{q), 
then the evolution equation for the Wigner function 
(Wigncr-Liouvillc equation) has the form 3i..l&|: 



dF^ 

dt 



m dq ^ ^ (2n-h 1)! ag2"+i a„2«+i 

71 — ^ ' ^ 



(10) 

If the potential does not have the terms with more than 
the second power of 5, then Eq. (|10|l has the same form 
as for a classical distribution function /: 

dJ_^p_dJ_^dV_dJ_ 
dt m dq dq dp 



(11) 
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The Wigner function must satisfy a number of conditions 
[Tsf . therefore initial function F^{q,p,t — 0) can not be 
chosen arbitrary. Even if F^{q,p,t) satisfies 'classical' 
equation (for specific potential V) it describes quan- 
tum system adequately because all quantum corrections 
(all powers of h) are held in the initial Wigner function 
F^{q,p,t ~ 0). For example, the uncertainty principle 
holds. 

One can rewrite (|10() in the form analogous to 



dF 



w 



p OF 



w 



dVeff 



dt m dq dq dp 

where a new effective potential 14/ / is introduced 

dVeff dF^ _ dV dF^ d^V d^F^ 



dp 



dq dp 24 dq-^ dp^ 



(12) 



(13) 



The characteristics of Eq. Hll|l obey the equations co- 
inciding with classical equations of motion 



dt 



P_, 

m' 



dp 



dV{q,P,t) 
dq 



(14) 



From Eq. (|12|) one can obtain the modified equations of 
motion for Wigner trajectories j3j 



dq 
di 



dp 
'dt 



9Veff{q,P,t) 
dq 



(15) 



To get information about time evolution of the system 
we numerically solve equations H15|) for the ensemble of 
trajectories. To simplify our calculation of Veff (|13() for 
the problem of interest we choose the analytical form of 
the external potential and the interaction between par- 
ticles to contain only the 2-nd and the 4-th powers of 
coordinates. In this case only the first two terms in the 
r.h.s. of Eq. (|l()|l are non zero. As a result the total force 
is a sum of the usual classical force and the 'quantum' 
force F^"°-"'* ^ which is infinite series in general, but in our 
case contains only one term: 



rpquant 



g3y Q2pW ^ 

dqidqidqk dpidpi F^ ' 



(16) 



where index k is the fc-th component of the force vector 
(there are N x d such components, N is the number of 



particles and d is spatial dimensionality), repetition of 
indexes indicates the summation. 

As one can note the 'quantum force' depends on the 
Wigner function, which is unknown. To overcome this 
problem we use a local approximation for the Wigner 
function in the vicinity of phase space point Xa by Gaus- 
sian 0: 



F'^{q,p,t) = F,] 



W -[{x-Xa{t))A^{t){x-Xa{t)) + b^it){x~Xait))] 

(17) 

where x = (^) is vector of all particle coordinates and 
momenta, matrix Aa (in our case of dimensionality 4x4) 
and vector ba (with dimensions 4x1) are obtained from 
the local moments of the ensemble of trajectories in the 
vicinity of point Xa- 

B. Consideration of exchange 



Exchange effects in this method can be in some cases 
considered simply by using special initial conditions. 
Consider the system with the wave function '^{x,t). One 
can obtain the Wigner function as 



F'^{q,p,t) 



(2 



^ld^e^^^f-^*iq + l,t)^iq-l,t) 



(18) 

If the system consists of either bosons or fermions, wave 
functions must be symmetrical or anti-symmetrical. If we 
regard the case when the Hamiltonian does not depend 
on spins of the particles, then we can consider only the co- 
ordinate part of wave function. Depending on the overall 
spin the coordinate part of wave function is either anti- 
symmetrical or symmetrical. For example, wave function 
of the following form is symmetrical (anti-symmetrical): 



1^(1,2)) 



\Ml)) 102(2)) ±|0i(2)) 10,(1)) 



(19) 



where (i) means the dependence on variables of the i-th 
particle. We use this wave function for the initial system 
state with of the form of a Gaussian wavepacket. As 
a result the Wigner function takes the form: 



P2) 



2(1 ± 


(</>i|</>2)P)(27r;i)2 












01 (qi - y) 02 



and can be rewritten as 



(92- f) 02 (91 -y) 



F^{qi,q2,pi,P2) = , I ,^ I , [Wi{qi,pi)W2iq2,P2) + 

2(1 ± I (0i|02)r) 

Wi{q2,P2) 1^2(91, Pi) ± Ui2{q2,P2) U2liqi,Pl) ± ?7l2(gi,Pl) U2l{q2;P2)] , 



(20) 



(21) 
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where 

Wk{q,p) 
and 



UkMp) = / dx e'^^l^ + |) -^.(9 - |) (23) 

In coordinate space the initial state (|19|l is described 
by wave function of the following form: 



4>k{x) 



{2nh) 
For this case, 



1 e,p(J^_^ + !PMli^_^](^2i) 



W,{,,p)^-exp^ ) (25) 



and the term with Ukj in Eq. I|21|) is proportional to 
exp (— A(a;io — 2^20)^)1 where is a positive constant. If 
CTi = 0-2 = then ^ = l/(2cr^). For |xio — a;2o| >> tr 
this term can be neglected and one gets: 



F^{qi,q2,Pi,P2) = 



1 



2(l±|(0i|02)P) 

[M^i(qi,pi) W2iq2,P2) + Wi{q2,P2) W2iqi,pi)] (26) 

We emphasize that this approximation is used only at the 
initial time moment. Further the dynamical equations 
are solved formally exactly. 

In the considered problem two particles move in the 
potential 



U{x) = «(— + 72:^), a, 7 > 



(27) 



The potential of interparticle interaction is: 
Vint = {const — Pr^}, if {const — Pr^} > 0, and = 0, oth- 
erwise. If we disregard discontinuity in the interparticle 
potential then the distinction from harmonic oscillator is 
the 4-th power of x and one has only one quantum term 
in the force pt)|) . Using the classical trajectories and 
the Gaussian approximation for the Wigner function one 
can solve the Wigner-Liouville equation 'exactly'. The 
distinction of the adopted approximation from the case 
of distinguishable particles is that now initial positions of 
two particles may be in the Gaussian centered at xiq or 
at X20 ■ In this way we regard the symmetry in exchange 
of particles and obtain the picture of their motion. 



C. Algorithm and calculation of average values 

Our simulation algorithm is the following. First, the 
initial coordinates and momenta of every trajectory in 
the ensemble are distributed according to chosen param- 
eters of the wavepackets (mean coordinate, momentum 
and their variances). Second, we calculate the 'quantum 
force' and solve numerically equations of motion. 



For the j-th trajectory at time t with coordinates and 
momenta {q^^\t), p^^\t)} one has to compute the lo- 
cal moments of the ensemble of trajectories in the vicin- 
ity of the point {q^^^t), p^^\t)} (point Xa), using the 
weight function which rapidly goes to zero with the in- 
crease of distance to this point (uncertainty principle 
must hold) 0. The approximation (|17|l is the many- 
dimensional Gauss distribution of the vector x = (^ . 

Matrix Aa — \C~^ , is the inverse matrix of covari- 
ance Ca, and vector ha = —2Aafa, where fa is the vector 
of averages. If one calculates (g^), (pi), {qtqk}, {PiPk), 
{qiPk), i,k ^ l,...Nd, one obtains fa = and 
symmetrical matrix Ca with elements 



Cail,m) = {{qi-ql''^){q, 



{qi-ql''^){q 



-qin- 



(28) 



for ; = 1,. ..iVd, m = 1,. 

Cail,m) ^ {{pi-Nd - p\%d)iPm-Nd - p';^lNd)) " 

{{Pl-Nd- p'f^Ndfl^iPrn-Nd- P^^^Nd)) (29) 

for / = A^d + 1, . . . 2Nd, TO = iVd -t- 1, . . . and 

Ca{l,m) = {{qi - q\"'^){Pm-Nd ~ P^^^Nd)) - 

{{qi~q\''^)){{p^-Nd-p^l^,)) (30) 

for I — 1, . . . Nd, TO = Nd+1, . . . 2Nd. Here qi andpi are 
the i-th components of vectors of all coordinates q and 
momenta p; (...) means the averaging over all trajecto- 
ries with the weight function, which rapidly approaches 
zero with growing distance to Xa in phase space. After 
that, one can calculate the inverse matrix for Ca and get 
Aa and ba- 

At every time t for the ^'-th trajectory with coordinates 
q^-^^t) and momenta p^^\t) matrix A^^^t) and vector 
b'^^\t) are calculated, therefore the 'quantum force' for 
the j-th trajectory is known. Further, one has to solve 
equations (|15|l . for example by Runge-Kutt and Adams 
methods. 

To calculate average values we use the following ap- 
proximation for the Wigner function: 

1 ^ 

F'^iq.p, t) = -Y.^{q- qk{t)} 5{p - pu{t)} (31) 



fc=l 



The summation is over all trajectories in the ensemble. 
One of the interesting values characterizing tunneling is 
the reaction probability: 



R{qa,t) 



1 



p{x, t) dx, 



(32) 



qa 



where the lower limit of integration is the point of the 
largest height of the barrier, p{x, t) is the particle density 
at point X and A'^ is the particle number. This quantity 
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shows what part of wavepackets are currently in the right 
weU. 

Another important value is f,unneling time. In this pa- 
per to determine tunneling time we use two methods: 
presence and arrival times. In the Wigner formalism 
these quantities are expressed through the following in- 
tegrals (see (Pi and (@J): 



|V'(a:o,t)|" = / dpF'^{xo,p,t) 



(33) 



and 



J{xo,t) 



J dppF'^ix„,p,t) 



(34) 



dt I dppF^{xa,p,t) 



(the Weyl symbol for flux operator has the form 
J.„ = |sin(M2|^) p{q-x,)). 



III. MODEL PROBLEM 

A. Hamiltonian 

We consider the following model problem. Two parti- 
cles move in one-dimensional space {e.g. in a quantum 
wire). The Hamiltonian of the system reads: 

^ = E + + ^'?')) + U - , (35) 

where qi,q2,pi, P2 are particle coordinates and momenta, 
U is interaction energy. We use the system of units with 
{h = m = a = 1}, Iq — fi'^'^ jimct)^!^ is the unit of 
length, Eq = h{a/my^^ is the unit of energy, and the 
unit of time is to = {m/ay^^. 

Initially particles are placed in the left well, their wave 
functions have the form of Gaussian wavepackets. Initial 
mean momenta and coordinates of wavepackets and their 
variance in momentum and coordinate spaces are chosen 
to make the transmission above the barrier as lower as 
possible and the overlapping of the wavepackets is negli- 
gible. Particles move in the direction of the barrier. This 
model roughly describes nonstationary tunneling of two 
electrons through the potential barrier in a quantum wire 
or tunneling between two quantum wells. This can be 
realized, for example, when with the help of laser pulses 
one prepares the state of two electrons in the form of 
two wave packets in nanostructure and study the system 
evolution in time. 

In the used system of units the Hamiltonian is: 

H-T.(4-l'+ + UiWi 92|) (36) 
i=i ^ ^ 



The Coulomb potential describes the inter- 

action between particles. The problem becomes one- 
dimensional if characteristic energies of the transverse 
quantization are much larger than the energies of the 
longitudinal motion. If it is valid then adiabatic approxi- 
mation applies and the problem is really one-dimensional. 
The interparticle interaction then reduces to 



^ 



-.dp ■ 



X 



^(l-er/(-)). 



(37) 



where we performed integration over the particle wave- 
functions of transverse quantization, with a being the 
characteristic width of the quantum wire. The interac- 
tion parameter, A = QxQ^m?^^ lij^^'^oi^l'^)^ is the ratio 
of characteristic Coulomb energy and energy of oscillator. 

In our model problem we substitute the potential H37|) 
for model quadratic one, just to avoid uncertainties re- 
lated with the calculation of the quantum force H16() and 
demonstrate the method for the case of exchange. Pa- 
rameters of this potential are chosen to make it as close 
as possible to expression (|S7|l : 



C/(r) 



A(V^-0.05(r/a)2), if r < a(2O0F)i/2, 
0, if r > a(20V^)i/2. 



(38) 



B. Initial parameters 



We analyze the system described by the Hamiltonian 
(|36|l for two cases - with and without exchange, respec- 
tively. The main quantity analyzed is the reaction proba- 
bility H32|l . Its largest value is unity when both particles 
are entirely in the right well. The reaction probability 
clearly shows the 'distribution' of particles between two 
wells and characterize their time evolution. 

Consider two cases, we call them 0-t/i order WMD 
(Wigner molecular dynamics in zeroth approximation) 
and n-th order WMD (further simply refered to as the 
^quantum case^). For 0-th order WMD the quantum term 
in force l|lt)|l is neglected and evolution of the trajectories 
is determined by classical Hamilton equations. There- 
fore only passing above the barrier is taken into account. 
This approximation is not purely classical, because the 
initial distribution is the same as for the ^quantum case^: 
\F^ {q,p,t = 0)1 can contain arbitrary powers of Ti. As a 
result we have a classical evolution of the quantum dis- 
tribution function, therefore we call this case 0-th order 
WMD, not the 'classical' MD. The difference between 
0-th order and n-th order is that in the latter case the 
quantum term in the force is regarded. 

The initial distribution in coordinate space has the 
form of two Gaussian wavepackets, which practically do 
not overlap. In the ^quantum case^ one can consider two 
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situations. First, we can consider the problem neglect- 
ing exchange (i.e. regard the distinguishable particles). 
Second, we can take exchange into account. In the first 
case the initial wave function is a product of onc-particle 
wave functions, and the Wigner function has the form 
{qi,q2,pi,p2) ^ Wi{qi,pi)W2{q2,P2) (compare with 
Eq. H2t)|) '). This means that one of the Gaussians corre- 
sponds to the first particle and another to the second one. 
For such initial distribution function exchange effects can 
not arise and we will call this situation the ' quantum case' 
without exchange. 

In the second situation the particles are identical and 
the initial wave function is symmetrical (or antisymmet- 
rical). Now the Wigner function has the form (|26|l . 
Both Gaussians may correspond either to the first or 
to the second particle. The initial coordinates and mo- 
menta of some trajectory, {xi,X2, 2/1,2/2}, are chosen with 
the probability {qi,q2,pi,p2)\ of the configuration 
{qi = X2, 92 = a^i,Pi = 2/2,P2 = 2/i}- If the wavepackets 
are initially close to each other, the terms Ukj in Eq. 1)21(1 
do not vanish and the procedure of setting the initial co- 
ordinates and momenta becomes more complicated. 

For both ^quantum cases' (without and with exchange) 
dynamical correlations are taken into account due to so- 
lution of the Wigner-Liouville equation. Statistical cor- 
relations are not regarded in the ^quantum case' with- 
out exchange but they are allowed in the case with ex- 
change. In this sense two situations resemble Hartree and 
Hartree-Fock approximations, respectively. Note that we 
do not use mean-field approximation, the similarity must 
be regarded only in the meaning formulated above. 

For the Q-th order WMD both ways of setting the ini- 
tial coordinates and momenta for trajectories can be ap- 
plied, but it was found that the result is practically in- 
dependent on it. We use the following parameters: ini- 
tial coordinates and momenta for Gaussians are xi (0) — 
-140, pi(0) = 45 and a;2(0) = -310, p2(0) = 90, dis- 
persions in coordinate space {ax for the Gauss function 
ea;p[— x^/(2ct^)]) are the same for both wavepackets and 
equal 20, in momentum space Cp = fi/{2(Jx) = 0.025. 
Parameters of the external potential are a = 1, 7 = 
1.25 X 10-8. 

One of our aims is to investigate the influence of in- 
teraction on tunneling. We change the interaction pa- 
rameter A, starting with A = (no interaction). Effec- 
tive charge of electrons or holes in nanostructures can 
be controlled by changing the permittivity, but not in 
very wide range. Here we would like to draw the atten- 
tion to the fact that due to the use of special system of 
units, h = m — a — 1^ the region of variation of the 
dimensionless interaction parameter can be pretty wide. 
Actually, in this unit system the parameters of external 
potential are used. Therefore we can vary the interaction 
A by changing the external potential. This change leads 
to scaling of the units of length, time, energy and so on. 



IV. RESULTS AND DISCUSSION 

The definition of tunneling is usually given for one par- 
ticle. We analyze the motion of two particles and as the 
under barrier transmission is of interest for us, the total 
initial energy of the system is set about 0.99 of the height 
of the barrier. So even if one of particles borrows all the 
energy the latter is still lower than the height of the bar- 
rier. We deal with the wavepackets, therefore though the 
mean energy per particle is lower than the barrier height, 
the transmission above the barrier is still possible due to 
the dispersion in the momentum space. The comparison 
between {)-th and n-th order WMD allow us to estimate 
the transmission under the barrier. The state of the par- 
ticles is considered as an entangled unified whole and the 
'transmission under the barrier' means the transmission 
of at least one of the particles. 



A. Reaction probabilities 

In Fig. ^ we show the reaction probability for the in- 
teraction A 0, 2 X 10"*, 6 X lO'' and 2 x 10^. The 'quan- 
tum case^ (with exchange) and 0-th order WMD are com- 
pared. One can see that for weak interaction there are 
a large difference between this two cases. Under barrier 
transmission takes place only for the 'quantum case', in 
the 0-th order only wavepacket components with the en- 
ergy above the height of the barrier can pass to the right 
well. 

With the increase of the interaction the reaction prob- 
ability grows for both cases. The reason is that the initial 
energy becomes larger with the increase of A and there are 
more wavepacket components with the energy above the 
barrier height. For the 'quantum case' it is also impor- 
tant that there are some high-energy components which 
pass under the barrier. For very large values A > 6 x 10* 
reaction probabilities of the 0-th order and the 'quantum 
case' are almost the same. It means that for strong in- 
teraction the role of tunneling is negligible, wavepackets 
include too many components, which can pass above the 
barrier. 

In Fig. |2] we present reaction probabilities for the 
'quantum case' with and without exchange. For large A 
classical transmission above the barrier prevails and the 
influence of exchange on tunneling is considerably small. 
For small A one can see a large difference in the reaction 
probabilities. If the particles are distinguishable, the re- 
action probability is larger, but the sign of the effect de- 
pends on the initial parameters of the wavepackets, they 
can be fitted to make reaction probability greater for the 
case with exchange. Initially the distribution of the par- 
ticles in coordinate space has the form of two separated 
Gaussians. We found that, for the used parameters, these 
two peaks quickly merge, forming a single peak, and very 
seldom they can be seen as two separated wavepackets. It 
means that particles are close to each other during most 
time of simulation and exchange effects must be substan- 
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FIG. 1: Time-dependence of the reaction probability, -R(t), for the Q-th order WMD (dashed lines) and ^quantum case! (solid 
lines). Interaction strengths are A = 0, 2 x 10'', 6 x 10'', 2 x lO''. Maximum value of the reaction probability R(i) = 1, time t 
is in the units to = {rajaY^'^ , a is the potential parameter l|35|l . 



tial. From Fig. |21one can see that exchange almost does 
not influence the reaction probability for large A. In this 
situation the contribution of transmission above the bar- 
rier to the reaction probability is very high (in compar- 
ison with tunneling, compare with the 0-th order WMD 
in Fig. So it is difficult to notice exchange effects 
against this background. It is possible that the observed 
effect of the increase of transmission above the barrier 
with the growth of the interaction parameter masks the 
effect of the increase of dynamical correlations between 
electrons, which also suppress exchange effects. 

The reaction probabilities presented in Figs. ^ and [3 
were measured with the finite precision, connected with 
the finite time step and finite number of the used trajec- 
tories. Therefore the regions of constant reaction prob- 
ability for small A arise - the subtle noise was smoothed 
over the errors. The substantial features (differences) of 
0-th order and n-th order WMD in cases with and with- 
out exchange become apparent on much larger scales and 
they are easily taken into account. Probably, the ob- 
served noise is due to the additional oscillations of high- 
energy parts of the wave packets. It will be the subject 
of detailed investigations in the next work. 

Another common feature of Figs.^E^ndl^lis the behav- 
ior of the reaction probability for A = 2 x 10^. At the end 
of considered time interval, < t < 30, the reaction prob- 
abilities settle around the value « 0.15. It seems that in 
the limit of infinite A the equilibrium corresponds to the 



situation when both the wells are occupied with the equal 
probability, because in this case the total initial energy 
is much greater than the height of the barrier. Then the 
reaction probability must settle around 0.5. Probably, 
A = 2 x 10^ is not large enough and the value of ini- 
tial total energy is sufficient only to push a small part of 
the wave packet through the barrier. It is possible that 
then the leakage to the right well is balanced by the op- 
posite transmission to the left well and that is why the 
reaction probability, R{t), stays near the value 0.15. On 
the other hand, there can be another explanation: after 
the transmition of some part of the wavepackets to the 
right well, the repulsion between this part and the part 
in the left well prevents further penetration of particles 
to the right well, as a result the reaction probability set- 
tles at « 0.15 (i. e., mechanism analogous to Coulomb 
blockade). Whether this second mechanism realizes is 
unclear, perhaps, these two processes take place simulta- 
neously. As for the first scenario, the kinetic equilibrium, 
there are no doubts it can exist: the reaction probability 
also settles around some value for A = (the ^quantum 
case' without exchange. Fig. [21). For this case there is 
no interaction, therefore the only explanation can be the 
equality of transmissions through the barrier in both di- 
rections. In fact, there are no reasons to believe that 
the reaction probability will always stay near, say, 0.15 
(for A = 2 X 10^), possibly here we observe only the in- 
termediate equilibrium and later the system can come 
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FIG. 2: Time-dependence of the reaction probability, R{t), for the 'quantum case' with exchange (solid lines) and without 
exchange (dashed lines). Interaction strengths are A = 0, 2 X 10*, 6 X lO", 2 X 10^ 



to the state when the reaction probabihty is about 0.5. 
But such extra-long-time evolution must be the subject 
of special investigation. 



B. Hartree description of the tunneling 

Both Fig. Hand Fig. [3 demonstrate oscillations of the 
reaction probability, with the growth of interaction pa- 
rameter A their period decreases and the picture becomes 
less regular. This is due to the behavior of the transmit- 
ted part of wavepackets: to the right from the barrier 
there is the wall of the right well, transmitted part is re- 
flected from it and moves to the left. Then it is partially 
transmitted back to the left well, making some modula- 
tion of the reaction probability curve. Transmission takes 
place mainly when the wavepackets come to the barrier, 
therefore reaction probability changes step by step with 
some period. When the interaction is weak, the particles 
move almost independently and this period coincides ap- 
proximately with that of oscillations of one particle in the 
well (T « 4, see Figs.n]and|3 the curves for A = 0). But 
with the increase of A this period is changed accordingly 
to the influence of one particle on the motion of the other 
particle. This influence can be illustrated more transpar- 
ently in the picture of effective one-particle barriers (see 
below). 

In terms of Hartree approximation it means that each 



particle moves effectively in the double-barrier potential. 
The first barrier is the stationary barrier between the 
wells and the second one is the effective time-dependent 
barrier due to the interaction of particles. Each particle 
either fall on two barriers, or move between them. In 
the latter case a particle is inside the potential well. As 
the other particle moves closely, the well becomes more 
narrow and the energy levels in it higher. Therefore the 
energy of the particle in the well becomes greater and 
the probability of transmission increases. Another mech- 
anism is that due to the nonadiabatic narrowing of the 
well the particle can jump on the higher energy levels, 
which also results in the increase of transmission proba- 
bility. 

The effective potential for the first particle in Hartree 
approximation can be defined as 

Veff{xi,t) = Vext{xi) + j U{xi,X2)\il{x2,t)\^ dX2, 

(39) 

where Vext is the external potential, U is the interaction 
between particles, and |'0(a;2,i)P is the probability den- 
sity for the second particle. The effective potential for 
the second electron is given by the analogous equation. 
Of course, this can be applied only to the case when the 
particles are distinguishable. If the particles are iden- 
tical, one can use the Hartree-Fock approximation, and 
the potential becomes nonlocal. We use neither Hartree 
nor Hartree-Fock approximation in our method, they are 
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FIG. 3: Effective potentials 1)39^ (solid lines), compared with the stationary external potential (dashed lines), in units of the 
height of the barrier. Two times, t = and t = 0.5, are considered. Two left (right) plots are the effective potentials for the 
particle, which is initially closer (farther) to (from) the barrier. We consider the 'quantum case' without exchange, A = 2 x 10^, 
X is in units Iq = h^^'^ (ma)"^^'^ , a is the potential parameter (I35II . 



just very convenient tools of visualizing the behavior of 
quantum particles, which interact with the barrier and 
between each other. 

In Fig. O we plot the effective potentials (|39|l . A = 
2 X 10^. The dotted line is the shape of external poten- 
tial. The 'quantum case' without exchange is considered, 
the particles are not identical, therefore Hartree approx- 
imation can be used. The case with exchange is more 
interesting but if particles are identical, one can regard 
the effective potential, which is the same for all parti- 
cles. Here we just illustrate the possible mechanisms of 
transmission and the case of distinguishable particles is 
more representative. In this case one can differ two situa- 
tions. First, the barrier, which arises due to interparticle 
interaction, is close to the the stationary barrier (effec- 
tive broadening of the stationary barrier in the two right 
plots in Fig. inj. Second, this effective barrier is closer 
to the left wall of the well (the two left plots in Fig. O . 
For identical particles, these situations take place simul- 
taneously and the effective potential is the same for every 
particle. 

The right plots in Fig. |2| show that the interaction 
makes the barrier wider and it prevents the transmis- 
sion. But in the left plots, the well becomes not so deep, 
the energy levels in it grow and low-energy components of 
the wavepacket can be reflected from the effective barrier 
in the direction of the stationary barrier. The tunnel- 



ing and transmission above the barrier from the higher 
energy levels are stronger. The oscillations of the low- 
energy components between the stationary and effective 
barriers make the tunneling probability larger. 

C. Tunneling times 

Consider now Table^where we present tunneling times 
for the spatial interval [—45, 45] determined by two meth- 
ods: presence ^ and arrival times (0J. The interaction 
strengths are A = 0, 6 x 10*, 2 x 10^. We consider two sit- 
uations: the 'quantum case' with and without exchange. 
At the edge points of this spatial interval the value of 
the external potential approximately coincides with the 
initial energy of the system. Two methods give close re- 
sults, tunneling time calculated with the help of Eq. Q 
is in general greater than that with use of Eq. but 
the trend connected with changes in tunneling time with 
change of the interaction parameter is the same for both 
methods. Tunneling time is greater for the case with 
exchange, because for our parameters there are fewer 
high-energy components in wavepackets and tunneling 
is weaker for this case. With the increase of A the tun- 
neling time gets smaller, this can also be connected with 
the fact that for strong interaction there are more high- 
energy components in the wavepackets. Those compo- 
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nents make the main contribution to the transmission 
because the more the energy the greater the transmis- 
sion probabiUty. Because the components with higher 
energy move faster the time of passing some interval be- 
come smaller. We do not list in Table transmission 
times for the case of free motion. These times are all 
about 2.0 and are almost independent on A, exchange 
and method of their calculation. One can see that the 
presence of the barrier makes the transmission (tunnel- 
ing) time smaller. It is due to enriching of the transmit- 
ted part of the wavepacket by high-energy components, 
which has greater probability to go through the barrier 
and move faster. Therefore on average transmitted part 
arrives to detector earlier than the whole wavepacket in 
the case of free motion. 



TABLE I: Tunneling times for the spatial interval [-45,45]. 
(The system of units is?i = m = a = l). 





A = 


A = 6 X 10"' 


A = 2 X 10^ 


presence time (no exchange) 
arrival time (no exchange) 


0.8(±0.1) 
0.9(±0.1) 


0.72(±0.09) 
0.86(±0.09) 


0.38(±0.05) 
0.80(±0.08) 


presence time (exchange) 
arrival time (exchange) 


0.8(±0.1) 
1.2(±0.1) 


0.78(±0.08) 
1.06(±0.09) 


0.36(±0.04) 
0.72(±0.07) 



ics method based on the Wigner representation. The 
WMD method allows to calculate different features char- 
acterizing quantum evolution and influence of particle 
interaction and exchange on tunneling. We found that 
the strong interaction in this problem leads to decrease 
of the role of tunneling in the transmission. If the inter- 
action is not very strong (A < 6 x 10^) then exchange 
effects are substantial and affect tunneling. The devel- 
oped WMD method allowed us to analyze some interest- 
ing features of tunneling and proved to be a powerful tool 
for study of nonstationary quantum processes. Of course, 
the numerical solution of Schrodinger equation would not 
be too difficult for the problem under consideration, but 
in this work we just demonstrate the method, regarding 
relatively simple system of two identical particles in the 
double-well potential. We made only the first step in the 
development of WMD for the many-body problems and 
we intend to report the results of investigation of many- 
particles system in the next paper. 



CONCLUSION 
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